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ABSTRACT 

We  investigate  the  effect  of  the  accelerated  expansion  of  the  Universe  due  to  a  cosmological 
constant,  A,  on  the  cosmic  star  formation  rate.  We  utilize  hydrodynamical  simulations  from 
the  eagle  suite,  comparing  a  ACDM  (cold  dark  matter)  Universe  to  an  Einstein-de  Sitter 
model  with  A  —  0.  Despite  the  differences  in  the  rate  of  growth  of  structure,  we  find  that  dark 
energy,  at  its  observed  value,  has  negligible  impact  on  star  formation  in  the  Universe.  We  study 
these  effects  beyond  the  present  day  by  allowing  the  simulations  to  run  forward  into  the  future 
(t  >  13.8  Gyr).  We  show  that  the  impact  of  A  becomes  significant  only  when  the  Universe 
has  already  produced  most  of  its  stellar  mass,  only  decreasing  the  total  comoving  density  of 
stars  ever  formed  by  15  percent.  We  develop  a  simple  analytic  model  for  the  cosmic  star 
formation  rate  that  captures  the  suppression  due  to  a  cosmological  constant.  The  main  reason 
for  the  similarity  between  the  models  is  that  feedback  from  accreting  black  holes  dramatically 
reduces  the  cosmic  star  formation  at  late  times.  Interestingly,  simulations  without  feedback 
from  accreting  black  holes  predict  an  upturn  in  the  cosmic  star  formation  rate  for  t  >  15  Gyr 
due  to  the  rejuvenation  of  massive  (>10n  Mq)  galaxies.  We  briefly  discuss  the  implication 
of  the  weak  dependence  of  the  cosmic  star  formation  on  A  in  the  context  of  the  anthropic 
principle. 
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1  INTRODUCTION 

Precise  observational  data  from  the  past  two  decades  have  allowed 
us  to  measure  the  cosmic  history  of  star  formation  back  to  very 
early  times  (z  ^  8).  The  star  formation  rate  (SFR)  density  of  the 
Universe  peaked  approximately  3.5  Gyr  after  the  Big  Bang  (z  &  2), 
and  declined  exponentially  thereafter  (for  a  review,  see  Madau  & 
Dickinson  2014). 

Galaxy  formation  and  evolution  are  a  highly  self-regulated  pro¬ 
cess,  in  which  galaxies  tend  to  evolve  towards  a  quasi-equilibrium 
state  where  the  gas  outflow  rate  balances  the  difference  between 
the  gas  inflow  rate  and  the  rate  at  which  gas  is  locked  up  in  stars 
and  black  holes  (BHs)  (e.g.  White  &  Frenk  1991;  Finlator  &  Dave 
2008;  Schaye  et  al.  2010;  Dave,  Finlator  &  Oppenheimer  2012). 
Consequently,  the  cosmic  SFR  density  is  thought  to  be  determined 
both  by  the  formation  and  growth  of  dark  matter  haloes  and  by  the 
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regulation  of  the  gas  content  in  these  haloes.  The  former  depends 
solely  on  cosmology,  whereas  the  latter  depends  on  baryonic  pro¬ 
cesses  such  as  radiative  cooling,  stellar  mass  loss,  and  feedback 
from  stars  and  accreting  BHs. 

Which  of  these  factors  is  most  responsible  for  the  decline  in 
cosmic  star  formation?  It  could  be  driven  by  the  ‘freeze  out’  of  the 
growth  of  large-scale  structure,  caused  by  the  onset  of  accelerating 
cosmic  expansion.  As  galaxies  are  driven  away  from  each  other  by 
the  repulsive  force  of  dark  energy,  accretion  and  merging  are  slowed 
and  galaxies  are  gradually  starved  of  the  raw  fuel  for  star  formation. 
Or,  it  could  be  caused  primarily  by  the  onset  of  efficient  stellar  and 
BH  feedback. 

The  discovery  of  the  accelerating  expansion  of  the  Universe  was 
a  breakthrough  achievement  for  modern  cosmology  (Riess  et  al. 
1998;  Perlmutter  et  al.  1999).  However,  the  driving  force  behind  the 
acceleration  (generically  known  as  dark  energy)  is  still  unknown. 
At  present,  all  cosmological  observations  are  consistent  with  a  cos¬ 
mological  constant,  or  a  form  of  energy  whose  density  remains 
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constant  as  the  Universe  expands.  One  such  form  of  energy  is  vac¬ 
uum  energy:  the  energy  of  a  quantum  field  in  its  ground  state  (zero 
particles).  The  present  best-fitting  cosmological  model,  known  as 
the  concordance  model,  or  ACDM,  includes  both  a  cosmological 
constant  A  and  cold  (i.e.  non-relativistic)  dark  matter.  This  model 
has  been  very  successful  in  matching  the  observational  data. 

Nevertheless,  the  model  raises  a  number  of  fundamental  prob¬ 
lems.  Predictions  from  quantum  field  theory  for  the  vacuum  energy 
density  overestimate  the  observed  value  of  A  by  many  orders  of 
magnitude  (for  a  review,  see  Weinberg  1989;  Carroll  2001).  In  ad¬ 
dition,  the  energy  density  of  matter  and  the  cosmological  constant 
are  within  a  factor  of  a  few  of  each  other  at  the  present  time,  making 
our  epoch  unusual  in  the  evolution  of  the  Universe.  This  is  known 
as  the  coincidence  problem.  These  problems  have  motivated  the 
search  for  alternative  models  of  dark  energy  and  modifications  of 
gravity  that  might  explain  the  acceleration  of  the  Universe  more 
naturally.  For  example,  quintessence  models  propose  that  the  den¬ 
sity  of  matter  and  dark  energy  track  each  other.  In  many  models, 
however,  fine-tuning  of  the  model  parameters  is  still  required  to 
explain  their  observed  similarity  (see  e.g.  Weinberg  2000). 

An  alternative  approach  is  therefore  to  explain  the  observed  value 
of  A  on  anthropic  grounds.  This  has  already  been  applied  very 
promisingly  to  the  coincidence  problem.  Since  the  coincidence  con¬ 
cerns  the  time  that  we  observe  the  universe,  the  nature  and  evolu¬ 
tion  of  observers  in  the  Universe  are  highly  relevant.  For  example, 
Lineweaver  &  Egan  (2007)  argue  that  the  production  of  planets 
in  our  Universe  peaks  when  matter  and  dark  energy  are  roughly 
coincident  (see,  however,  Loeb,  Batista  &  Sloan  2016). 

For  the  cosmological  constant  and  other  fundamental  parameters, 
anthropic  reasoning  requires  a  multiverse.  Many  models  of  infla¬ 
tion,  such  as  eternal  inflation,  imply  that  the  Universe  as  a  whole  is 
composed  of  a  vast  number  of  inflationary  patches  or  sub-universes. 
Each  sub-universe  inherits  a  somewhat  random  set  of  physical  con¬ 
stants  and  cosmic  parameters  from  a  wide  range  of  possible  val¬ 
ues.  Sub-universes  in  which  the  cosmological  constant  is  large  and 
positive  will  expand  so  rapidly  that  gravitational  structures,  such  as 
galaxies,  are  unable  to  form  (e.g.  Weinberg  1987;  Martel,  Shapiro  & 
Weinberg  1998;  Efstathiou  1995;  Sudoh  et  al.  2017).  Large  negative 
values  will  cause  the  universe  to  collapse  rapidly,  also  preventing 
the  formation  of  galaxies.  Only  sufficiently  small  values  of  A  will 
lead  to  the  formation  of  universes  that  are  able  to  host  observers. 
This  argument  eliminates  extreme  values  of  A.  For  example,  Wein¬ 
berg  (2000)  estimates  an  upper  bound  on  a  positive  vacuum  energy 
density  to  allow  for  the  formation  of  galaxies  of  about  200  times 
the  present  mass  density. 

Refining  Weinberg’s  estimate  requires  us  to  more  accurately  ex¬ 
plore  the  sensitivity  of  galaxy  formation  to  the  presence  of  A.  Here, 
we  use  a  suite  of  hydrodynamical  simulations  to  take  a  first  look 
at  this  problem  by  calculating  the  effect  of  the  cosmological  con¬ 
stant  on  galaxy  and  star  formation  in  our  Universe.  Specifically,  we 
compare  the  formation  of  galaxies  in  our  Universe  with  a  hypothet¬ 
ical  universe  that  is  indistinguishable  from  ours  at  early  times  but 
has  no  cosmological  constant.  Because  A  is  negligibly  small  in  the 
early  universe,  these  two  universes  will  evolve  in  nearly  identical 
ways  for  the  first  r«2  Gyr  of  cosmic  time  (when  the  dark  energy 
density  is  less  than  0.03  times  the  matter  density).  This  means  that 


the  epochs  of  nucleosynthesis,  recombination,1  and  reionization  are 
indistinguishable. 

In  recent  years,  the  accuracy  of  our  understanding  of  galaxy  for¬ 
mation  has  improved  considerably,  reaching  the  point  at  which  it  is 
possible  to  undertake  this  comparison  meaningfully.  The  increased 
realism  of  simulated  galaxies  (in  particular  disc  galaxies  with  more 
realistic  sizes  and  masses)  has  been  achieved  due  to  the  use  of 
physically  motivated  subgrid  models  for  feedback  processes  (e.g. 
Schaye  et  al.  2015;  Dubois  et  al.  2016;  Pillepich  et  al.  2017).  One  of 
the  key  ingredients  that  have  allowed  this  progress  is  the  inclusion 
of  realistic  models  for  the  impact  of  feedback  from  the  growth  of 
super  massive  BHs  (e.g.  Bower  et  al.  2017).  All  successful  models 
now  demonstrate  the  need  for  active  galactic  nuclei  (AGNs)  as  an 
additional  source  of  feedback  that  suppresses  the  formation  of  stars 
in  high-mass  haloes  (e.g.  Benson  et  al.  2003;  Bower  et  al.  2006; 
Croton  et  al.  2006;  Crain  et  al.  2015;  Pillepich  et  al.  2017).  One  of 
the  aims  of  this  paper  is  to  compare  the  impact  of  the  cosmological 
constant  with  that  resulting  from  the  inclusion  of  BHs  in  the  simu¬ 
lation.  In  a  previous  study,  van  de  Voort  et  al.  (2011)  found  that  by 
preventing  gas  from  accreting  on  to  the  central  galaxies  in  massive 
haloes,  outflows  driven  by  AGN  play  a  crucial  role  in  the  decline  of 
the  cosmic  SFR. 

Different  groups  have  used  hydrodynamical  simulations  to  study 
the  effect  of  different  dark  energy  or  modified  gravity  models  on  cos¬ 
mological,  galactic,  and  sub-galactic  scales  (e.g.  Puchwein,  Baldi 
&  Springel  2013;  Penzo  et  al.  2014,  2016).  Taking  a  different  ap¬ 
proach,  in  this  paper  we  investigate  the  effect  of  the  accelerated 
expansion  of  the  Universe  on  galaxy  formation  by  asking  the  fol¬ 
lowing  question: 

How  different  would  the  Universe  be  if  there  had  been  no  dark 

energy? 

For  our  study,  we  use  a  suite  of  large  hydrodynamical  simu¬ 
lations  from  the  Evolution  and  Assembly  of  GaLaxies  and  their 
Environment  (eagle)  project  (Crain  et  al.  2015;  Schaye  et  al.  2015). 
Using  state-of-the-art  subgrid  models  for  radiative  cooling,  star  for¬ 
mation,  stellar  mass  loss,  and  feedback  from  stars  and  accreting 
BHs,  the  simulations  have  reproduced  many  properties  of  the  ob¬ 
served  galaxy  population  and  the  intergalactic  medium  both  at  the 
present  day  and  at  earlier  epochs  (e.g.  Furlong  et  al.  2015,  2017; 
Lagos  et  al.  2015;  Rahmati  et  al.  2015,  2016;  Schaller  et  al.  2015; 
Trayford  et  al.  2015;  Bahe  et  al.  2016;  Rosas-Guevara  et  al.  2016; 
Segers  et  al.  2016).  Given  that  the  physics  of  the  real  Universe  is 
reasonably  well  captured  by  the  phenomenological  sub-grid  mod¬ 
els  implemented  in  the  simulations,  with  the  use  of  appropriate 
assumptions,  we  can  run  the  simulations  beyond  the  present  time, 
and  explore  the  consequences  of  our  models  for  the  future.  Fur¬ 
thermore,  the  simulation  re-scaling  strategy  developed  here  will 
be  used  in  a  companion  paper  (Barnes  et  al.  2018),  considering  a 
wider  range  of  A  values  and  determining  the  likelihood  distribution 
of  possible  A  values  conditioning  the  existence  of  observers. 

The  layout  of  this  paper  is  as  follows.  In  Section  2,  we  develop 
a  simple  analytic  model  of  the  cosmic  SFR  that  captures  the  sup¬ 
pression  due  to  a  cosmological  constant.  In  Section  3,  we  briefly 
describe  the  simulations  from  which  we  derive  our  results  and  dis¬ 
cuss  our  criteria  for  halo  and  galaxy  definitions.  In  Section  3,  we 
also  describe  our  motivations  to  run  our  cosmological  simulations 


'Of  course,  an  observer  in  a  A  =  0  universe  would  measure  a  different 
angular  power  spectrum  in  the  cosmic  microwave  background  after  13.8  Gyr, 
because  of  the  very  different  expansion  history  of  the  Universe  at  later  times. 
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into  the  future,  and  our  assumptions  in  doing  so.  Section  4  provides 
a  detailed  discussion  of  our  re-scaling  strategy  for  the  alternative 
cosmological  models.  In  Section  5,  we  explore  the  dependence  of 
the  star  formation  history  of  the  universe  on  the  existence  of  a  cos¬ 
mological  constant  and  the  presence  of  BHs.  We  also  explore  their 
impact  on  other  galaxy  population  properties,  both  up  to  the  present 
time,  and  into  the  future.  Finally,  we  summarize  and  discuss  our 
results  in  Section  6. 


2  A  SIMPLE  ANALYTIC  MODEL  FOR  THE 
COSMIC  SFR  DENSITY 


2.1  Comparing  different  cosmological  models 


The  star  formation  history  of  the  Universe  is  determined  by  the 
interplay  of  cosmic  expansion  and  the  time-scale  at  which  cold 
gas  can  turn  into  stars.  These  processes  happen  on  time-scales  that 
differ  by  several  orders  of  magnitude,  but  are  coupled  through  the 
accretion  rate  of  gas  on  to  gravitationally  bound  haloes.  The  aim 
of  our  paper  is  to  compare  theoretical  universes  in  which  the  star 
formation  time-scales  are  the  same,  but  the  cosmological  time- 
scales  vary.  We  need,  therefore,  to  be  careful  when  comparing  the 
different  models,  since  the  choice  of  coordinates  that  vary  with 
cosmological  parameters  will  obscure  the  similarities  of  the  models. 
In  particular,  the  expansion  factor  at  the  present  day,  ag,  is  often 
treated  as  an  arbitrary  positive  number,  and  it  is  common  practice 
to  set  ag  =  1.  In  this  paper,  we  need  to  take  a  different  approach 
since  we  want  to  compare  the  properties  of  different  universes  at  the 
same  cosmic  time  (measured  in  seconds  or  a  multiple  of  key  atomic 
transitions).  Assuming  a  common  inflationary  origin,  normalizing 
out  flo  is  not  appropriate,  since  the  expansion  factor  at  the  present 
day  (to  =  13.8  Gyr)  would  be  different  for  each  universe. 

We  still  need  to  define  a  scale  on  which  to  measure  the  size 
of  the  universes  we  consider.  Using  a  hat  notation  Q  to  denote 
quantities  in  our  observable  Universe,  we  set  a0  =  a(t0 )  =  1.  We 
want  to  emphasize  that  the  cosmological  models  that  we  consider  all 
start  from  very  similar  initial  conditions.  It  therefore  makes  sense 
to  normalize  them  to  the  same  value  of  the  expansion  factor  at 
an  early  time,  t\.  We  therefore  set  a\  =  a(t\)  —  d(fi).  We  choose 
dj  =  1/(1  +  127),  corresponding  to  a  redshift  of  2  =  127  for  a 
present-day  observer  in  our  Universe.2  At  this  moment,  the  age  of 
the  universe  is  t\  =  1 1.98  Myr.  This  applies  to  all  of  the  universes 
we  consider  since  the  cosmological  constant  term  has  negligible 
impact  on  the  expansion  rate  at  such  early  times. 

Although  time  (in  seconds)  is  the  fundamental  coordinate  that 
we  use  to  compare  universes,  it  is  sometimes  useful,  for  example 
when  comparing  to  observational  data,  to  express  time  in  terms  of 
the  redshifts  measured  by  a  present-day  observer  in  our  Universe, 
z.  We  convert  between  cosmic  time  t  (which  is  equivalent  between 
universes)  and  a  by  inverting  the  time-redshift  relation  for  our  Uni¬ 
verse: 


ap 

a(t) 


(1) 


It  is  important  to  note  that  z  is  not  the  redshift  that  would  be  mea¬ 
sured  by  an  observer  in  an  alternative  universe. 

In  this  paper,  we  will  focus  our  comparison  on  two  cosmological 
models,  a  standard  ACDM  universe  as  inferred  by  Planck  Col¬ 
laboration  XVI  (2014),  and  an  Einstein  deSitter  (EdS)  universe. 


2Z  =  127  was  the  reference  simulation's  starting  redshift. 


Table  1.  The  cosmological  parameters  for  the  eagle  simulations  used  in  this 
study.  ACDM  model  refers  to  parameters  inferred  by  Planck  Collaboration 
(2014).  EdS  refers  to  an  Einstein-de  Sitter  universe.  £2m,  £2a.  Oi>  are  the 
average  densities  of  matter,  dark  energy,  and  baryonic  matter,  respectively, 
in  units  of  the  critical  density  at  redshift  zero;  Hg  is  the  Hubble  constant, 
o%(t\ )  is  the  square  root  of  the  linear  variance  of  the  matter  distribution  at  the 
initial  cosmic  time  of  the  simulations  {t\  =  1 1.98  Myr)  when  smoothed  with 
a  top-hat  filter  of  radius  11.8  cMpc  (8/t-1  cMpc  for  a  ACDM  model),  ns  is 
the  scalar  power-law  index  of  the  power  spectrum  of  primordial  adiabatic 
perturbations,  and  Y  is  the  primordial  abundance  of  helium.  Values  in  bold 
show  differences  with  respect  to  the  ACDM  values. 


Cosmological 

Parameter 

ACDM  (Ref) 

EdS 

0.307 

1 

Da 

0.693 

0 

Db 

0.048  25 

0.157  17 

h  = 

0.6777 

0.3754 

//0/(100kms-1  Mpc-1) 
rrsRl) 

0.0083 

0.0083 

ns 

0.9611 

0.9611 

Y 

0.248 

0.248 

Assuming  both  cosmological  models  have  a  common  inflationary 
origin,  the  models  can  be  normalized  as  follows: 

(i)  For  the  ACDM  model  (see  Table  1),  we  set  ap  =  a(tp)  = 
1,  where  tg  —  13.82  Gyr  is  the  present-day  age  of  the  universe. 
At  time  t\  =  11.98  Myr,  a\  =  a(t\)  =  0.007813.  At  this  time,  the 
expansion  rate,  as  measured  by  the  Hubble  parameter  is,  H\  = 
H(h)  =  54,  377  kms-1  Mpc-1  =55.6Gyr-1. 

(ii)  We  require  the  EdS  model  to  have  the  same  early  expan¬ 
sion  history,  i.e.  a{t\)  =  0.007813  and  H(ti)  =  55.6Gyr-1.  In 
this  universe,  at  the  present  day  (i.e.  t  —  tg  —  13.82  Gyr)  the 
universe  has  a  size,  ap  —  a(tg)  =  0.8589  and  an  expansion  rate, 
H(t0)  =  0.0482  Gyr-1  =  47.16kms-1  Mpc-1. 

Fig.  1  shows  the  cosmic  scale  factor  as  a  function  of  time  for  the 
two  cosmological  models.  As  expected,  at  t  =  t0,  an  EdS  universe 
is  smaller  in  size  at  the  present  day,  as  the  cosmic  expansion  has 
not  been  accelerated  by  the  effect  of  A.  As  the  two  universes  evolve 
into  the  future,  the  size  differences  and  relative  expansion  rates 
grow,  e.g.  at  t  =  20  Gyr,  the  scale  factor  for  the  ACDM  models 
is  Ri25  per  cent  larger  than  for  the  EdS,  and  the  expansion  rate  is 
Ri50  per  cent  larger  for  our  Universe. 


2.2  Cosmological  expansion  history  as  a  function  of  time 


In  the  standard  model  of  cosmology  for  a  homogeneous  and 
isotropic  universe,  the  geometry  of  space-time  is  determined  by 
the  matter-energy  content  of  the  universe  through  the  Einstein  field 
equations  as  described  by  the  Friedmann-Lemaltre-Robertson- 
Walker  metric  in  terms  of  the  scale  factor  a(t)  and  the  curvature 
K,  yielding  the  well-known  Friedmann  equation. 


,  8  nG  Kc 1  Ac2 

-  =H\t)=—  p-  —  +  — 

\a )  3  a1  3 


(2) 


where  H{t)  is  the  Hubble  parameter.  As  the  inflationary  models 
predict  that  the  Universe  should  be  spatially  flat,  we  only  consider 
universes  with  no  spatial  curvature,  i.e.  K  =  0. 

The  density  of  equation  (2)  includes  the  contribution  of  non- 
relativistic  matter  and  radiation  (pm  and  pr).  The  radiation  content 
of  the  Universe  dominated  its  global  dynamics  at  very  early  times 
{a  — >■  0),  but  its  contribution  is  negligible  thereafter.  Ignoring  pr  and 
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Figure  1.  Cosmic  scale  factor  as  a  function  of  time  for  two  cosmological 
models.  The  model  for  the  cosmological  parameters  for  a  standard  ACDM 
universe  as  inferred  by  Planck  Collaboration  XVI  (2014)  is  shown  in  blue. 
An  Einstein-de  Sitter  universe  is  shown  in  orange.  Note  that  by  construction 
the  scale  factors  are  indistinguishable  when  the  universes  are  less  than  1  Gyr 
old.  The  power  series  approximation  of  equation  (10)  is  shown  with  a  dashed 
green  line. 


using  the  energy  density  at  an  arbitrary  time  t\ ,  equation  (2)  can  be 
written  as 


8ttG 


An,l  ( 


+ 


Ac1 


(3) 


where  pm>  i  is  the  matter  density  of  the  universe  at  t  =  t\,  and 
a\  =  a(ti).  We  choose  q  such  that  it  corresponds  to  a  suffi¬ 
ciently  early  epoch,  when  the  contribution  of  the  cosmological  con¬ 
stant  term  is  negligible.  As  discussed  in  the  previous  section,  at 
this  time  any  universe  closely  approximates  an  EdS  universe  and 
we  can  assume  that  a\  =  a\  and  pmi  =  pml  =>  pm,o(4o/'Zi)3  = 
Pm.oGo/ai)3-  Then,  equation  (3)  can  be  written  as 


2 


(4) 


Note  that  in  equation  (4),  the  evolution  of  the  scale  factor  for 
any  arbitrary  cosmology  is  written  in  terms  of  the  matter  density  of 
our  Universe  at  the  present  time  pm,o.  We  have  left  the  factor  of  <3o 
explicit  in  the  equation,  but  it  can  be  set  to  ao  =  1 ,  noting  that  flo  7^ 
1  for  any  cosmological  model  different  to  our  Universe. 

The  LHS  of  equation  (4)  has  units  of  time-2  and  we  will  later 
find  it  useful  to  represent  the  RHS  as  the  sum  of  two  time-scales. 
The  cosmological  constant  is  often  written  as  an  energy  component 
with  energy  density  pa  =  Ac2/8nG ,  however,  we  can  express  this 
as  a  time-scale  as  follows: 


Similarly,  the  matter  content  of  the  Universe  can  be  expressed  as  a 
time-scale, 


tm  = 


3 

87rGpo 


(6) 


Using  the  cosmological  parameters  for  our  Universe,  tA  =  ?a  = 
17.33  Gyr  and  tm  =  tm  =  26.04  Gyr.  For  an  EdS  universe,  tA  — »■  oo. 
Using  this  notation,  equation  (2)  can  be  written  as 


-  I  =  L, 

a  . 


~  I  +  tA\ 


(7) 


which  can  be  solved  analytically  to  express  the  expansion  factor  as 
a  function  of  time  and  the  parameters  tm  and  tA : 


ait)  = 


1 



2/3 


In  the  limit  tA  — >  oo  this  reduces  to  the  familiar  EdS  solution, 

-|  2/3 


lim  a(t)  = 

tA-fOO 


3  t 

2  t„ 


(8) 


(9) 


In  order  to  explore  the  significance  of  the  t/tA  term  more  clearly, 
we  can  expand  equation  (8)  as  a  Taylor  series: 


a(t)  i 


3  t 
2  tm 


2/3 


1  + 


1 


1 

80  \tA 


+  • 


(10) 


The  coefficients  of  the  series  decreased  rapidly  so  that  the  first  three 
terms  provide  a  good  approximation  up  to  t  =  2 tA  and  beyond. 
Fig.  1  shows  how  well  this  power  series  approximation  works. 


2.3  The  growth  of  density  perturbations 

In  the  standard  model  of  cosmology,  structures  such  as  galaxies  and 
clusters  of  galaxies  are  assumed  to  have  grown  from  small  initial 
density  perturbations.  Expressing  the  density,  p,  in  terms  of  the 
density  perturbation  contrast  against  a  density  background, 

p(x,  0  =  p(t)[l  +  «(x,  f>],  (ID 

the  differential  equation  that  governs  the  time  dependence  of  the 
growth  of  linear  perturbations  in  a  pressureless  fluid,  such  as  dark 
matter,  can  be  written  as  (for  a  review,  see  Peebles  1980;  Mo,  van 
den  Bosch  &  White  2010), 

d2S  a  d<5 

r+2-—  -4nGp8  =  0.  (12) 

d t2  a  dt 

The  growing  mode  of  equation  (12)  can  be  written  as 
8(t)  =  D(t)8(t0),  (13) 


where  D(t)  is  the  linear  growth  factor,  which  determines  the  nor¬ 
malization  of  the  linear  matter  power  spectrum  relative  to  the  initial 
density  perturbation  power  spectrum,  and  is  computed  by  the  inte¬ 
gral 


a 

D(t)  oc  — 
a 


f 


At' 

a2{t') ' 


(14) 


Using  the  hat  notation  as  before,  we  normalize  D(t)  so  that 


(i)  £>+(f0)  =  1 

(ii)  D(tl)  =  £»(fl) 


In  general,  the  growing  mode  can  be  obtained  from  equation  ( 14) 
numerically.  Fig.  2  shows  the  growth  factor  as  a  function  of  cosmic 
time  for  the  two  cosmological  models.  As  expected,  the  figure  shows 
that  linear  perturbations  grow  faster  in  an  EdS  universe  compared 
to  those  in  a  ACDM  universe. 

It  is  possible  to  gain  more  insight  by  integrating  the  power-series 
approximation  for  a{t)  from  equation  (10).  Expanding  the  solution 
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Figure  2.  The  linear  growth  factor  for  the  ACDM  and  EdS  cosmological 
models.  The  rates  are  all  normalized  such  that  D(to )  =  1,  for  the  ACDM 
model,  and  D(t\)  =  D{t\).  The  bottom  panel  shows  the  growth  factor  at 
a  given  time,  divided  by  the  growth  factor  for  ACDM.  The  presence  of 
a  cosmological  constant  suppresses  the  growth  of  structure  in  the  ACDM 
model  (blue)  compared  to  that  in  the  EdS  model  (orange).  The  power  series 
approximation  of  equation  (15)  is  shown  with  a  dashed  green  line. 


Figure  3.  The  relative  rate  of  growth  of  density  perturbations,  ^  ^  for  the 
ACDM  and  EdS  cosmological  models.  The  bottom  panel  shows  the  ratio, 
at  a  given  time.  The  presence  of  a  cosmological  constant  slows  down  the 
growth  of  structure  in  the  ACDM  model  (blue)  compared  to  that  in  the  EdS 
model  (orange).  The  power  series  approximation  of  equation  (16)  is  shown 
with  a  dashed  green  line. 

2.4  Impact  on  halo  accretion  rates 


again  as  a  power  series  in  ( t/t a),  retaining  the  leading  terms,  yields, 


D(t)  = 


3  t 
2  tm 


2/3 


dlKD 


1-0.1591 


+  0.0366 


(15) 


The  growth  rates  of  linear  perturbations  do  not  directly  predict 
the  growth  rates  of  haloes;  however,  we  can  directly  connect  the 
two  through  the  approach  developed  by  Press  &  Schechter  (Press 
&  Schechter  1974;  Bond  et  al.  1991;  Bower  1991;  Lacey  &  Cole 
1993).  Correa  et  al.  (2015)  showed  that  the  accretion  rates  of  haloes 
can  be  written  as  (see  also  Neistein,  van  den  Bosch  &  Dekel  2006), 


where  Ku  is  a  normalization  constant.  Requiring  £>(fo)  =  1  gives 
Kd  =  4.70  x  10~3  Gyr  2.  Fig.  2  shows  that  equation  (15)  provides 
a  good  approximation  up  to  t  =  tA . 

This  demonstrates  that  although  the  tA  term  slows  down  the 
growth  of  perturbations,  its  effect  is  less  than  10  per  cent  until  t  ~ 
tA  (0.1/0. 159 1)1//2  ~  0.8/a  corresponding  to  »*13.8  Gyr  (~r0)  in 
our  Universe. 

As  we  discuss  in  the  following  section,  the  quantity  of  funda¬ 
mental  interest  for  the  accretion  rate  of  dark  matter  haloes  is  the 
relative  rate  of  growth  of  density  perturbations,  ^ .  We  show  this 
for  the  numerical  solution  in  Fig.  3.  We  can  also  compute  the  rela¬ 
tive  growth  rate  by  differentiating  the  power-series  approximation 
of  equation  (15).  Retaining  the  lowest  order  terms,  we  find 


1  AD 
D dT 


2 
3 1 


^1  -  0.4773 


(1)4°., «5 


(16) 


This  expression  does  not  depend  on  the  constants  tm  or  KD  because 
we  are  focusing  on  the  relative  change  in  the  growth  factor.  The 
impact  of  the  cosmological  constant  term  is  relatively  large,  creating 
an  r«50  per  cent  increase  in  growth  rate  for  the  EdS  model  compared 
to  ACDM  when  t  R»  tA. 


1  &Mh  =  [2  ( K/D )  1  AD 

Mh  At  V  n  S(Mhyr-  ( qy  -  1)1/2  D  At  ’ 

where  Mh  is  the  halo  mass  and  S(Mh )  is  the  variance  of  the  den¬ 
sity  field  on  the  length  scale  corresponding  the  halo  mass.  is  a 
parameter  that  represents  a  threshold  in  the  linearly  extrapolated 
density  field  for  halo  collapse.  The  parameters,  q  and  y,  are  related 
to  the  shape  of  the  power  spectrum  around  the  halo  mass  M/,.  Ap¬ 
proximating  the  scale  dependence  of  the  density  field  as  a  power 
law,  S  =  S0Mhy,  Correa  et  al.  (2015)  find  S  ^  3.98,  y  ^  0.3,  and 
q  ^  3.16,  giving  [S(Mh)(qy  -  l>r1/2  ss  0.78  for  1012Mq  haloes. 
These  values  depend  only  on  the  initial  power  spectrum  (which  we 
assume  to  be  the  same  in  all  the  universes  we  consider)  and  do 
not  depend  on  the  cosmological  parameters.  This  formulation  thus 
neatly  separates  the  contribution  of  the  power-spectrum  shape  from 
the  cosmological  parameters.  We  are  therefore  able  to  assume  that  q 
and  y  are  the  same  for  all  the  universes  that  we  consider,  and  focus 
on  the  dependence  on  D(t). 

For  the  numerical  values  of  the  power-spectrum  parameters 
around  a  halo  mass  of  1012  Mq,  equation  (17)  reduces  to 


1  AMh  1  dD 

- -  =  1.0456—  — . 

Mh  At  D2  At 


(18) 
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Figure  4.  The  specific  accretion  rate  of  haloes  of  mass  M,,  =  I012  Mq, 
for  the  ACDM  and  EdS  cosmological  models.  The  bottom  panel 
shows  the  ratio  at  a  given  time.  The  presence  of  a  cosmological  constant 
slows  down  the  specific  accretion  rates  of  haloes  in  the  ACDM  model  (blue) 
compared  to  that  in  the  EdS  model  (orange).  The  power  series  approximation 
of  equation  (19)  is  shown  with  a  dashed  green  line. 


This  dependence  can  be  understood  as  the  combination  of  two  fac¬ 
tors.  The  first  reflects  the  relative  growth  rate  of  density  fluctuations 
^  ^ .  The  second  factor  of  1  ID  comes  from  the  rarity  of  haloes, 
reflecting  the  higher  growth  rate  of  fluctuations  in  the  tail  of  the 
density  field  distribution. 

Further  insight  can  be  gained  by  using  the  series  approximation. 
This  gives 


1  d  Mh 

Jf~dT 


566.61 

■JStWtH3 


This  explicitly  shows  how  the  presence  of  a  cosmological  constant 
modulates  the  halo  growth  rate.  In  our  Universe,  the  impact  of 
the  cosmological  constant  term  is  relatively  modest,  however;  at 
t  =  13.8  Gyr,  we  expect  the  difference  to  be  20  per  cent,  growing  to 
40  per  cent  at  t  =  20  Gyr. 

As  an  example,  in  Fig.  4  we  show  the  accretion  rate  of  haloes  of 
Mi,  —  1012  Mq,  both  numerically  and  using  equation  (19). 


2.5  Impact  on  the  SFR  of  the  Universe 


In  order  to  link  the  SFR  of  haloes  of  mass  Mh  to  their  accretion  rate, 
as  a  first  approximation,  we  assume  a  time-independent  galaxy 
specific  SFR  to  host  halo  specific  mass  accretion  rate  relation  (e.g. 
Behroozi,  Wechsler  &  Conroy  2013a;  Tacchella,  Trenti  &  Carollo 
2013;  Rodrlguez-Puebla  et  al.  2016), 


3  log  M* 
Mh/Mh  ~  Slog Mh 


(20) 


where  the  star  formation  efficiency  e,  of  haloes  of  mass  Mh,  is  the 
slope  of  the  stellar-halo  mass  relation.  From  this  equation,  the  star 
formation  as  a  function  of  halo  mass  can  be  written  as 


M*(M„)  =  e,(Mh)Mh, 


(21) 


where  et(Mh)  :=  e{Mh)  x  is  completely  defined  by  the 

stellar-halo  mass  relation.  As  there  is  no  a  priori  knowledge  of  the 
functional  form  of  et(M/,),  we  use  the  abundance  matching  results 
from  Behroozi,  Wechsler  &  Conroy  (2013b)  to  estimate  e*{Mh). 
The  efficiency  e*(M,,)  peaks  at  masses  similar  to  Milky  Way-sized 
haloes  (~1012  Mq)  and  falls  steeply  for  higher  and  lower  masses. 
e.t(Mh)  can  be  well  approximated  by  a  broken  power  law  as 


e*{Mh)  oc  < 


^  1012Mq 


Mh 


Mh 


1012M, 


o 


if  Mh  <  1012Mq 


if  Mh  >  1012Mq 


(22) 


At  low  masses,  SFR  is  suppressed  because  of  the  efficiency  of 
feedback  from  star  formation,  at  higher  masses  the  cooling  of  the 
inflowing  gas  is  suppressed  by  heating  from  BHs  (White  &  Frenk 
1991;  Benson  et  al.  2003;  Bower  et  al.  2006;  Haas  et  al.  2013;  Crain 
et  al.  2015;  Dubois  et  al.  2016;  Bower  et  al.  2017). 

In  order  to  complete  the  analysis,  we  need  to  combine  the  specific 
halo  mass  accretion  rate  with  an  estimate  of  the  halo  abundance. 

In  the  Press  &  Schechter  analysis,  the  comoving  abundance  of 
haloes  of  mass  M/,  at  time  t  is  given  by  (Press  &  Schechter  1974; 
Bond  et  al.  1991;  Bower  1991;  Lacey  &  Cole  1993), 


d  n(Mh,t)  p0  Scy  1  / _ S_f_\ 

d Mh  Ml  fTnS'O-  D  6XP  V  2 SD2  ) 


(23) 


where  we  have  assumed  that  the  density  power  spectrum  is  a  power 
law  with  exponent  y  and  written  the  comoving  density  of  the  Uni¬ 
verse  as  po  following  our  convention.  Note  that  we  compute  co¬ 
moving  densities.  At  the  same  cosmic  time,  the  different  expansion 
rates  will  result  in  different  physical  (proper)  halo  and  SFR  den¬ 
sities,  simply  because  of  the  more  rapid  expansion  of  the  ACDM 
cosmology. 

The  total  cosmic  SFR  density  is  given  by  the  integral  of  all  star 
formation  in  all  haloes, 


f  .  d n(Mf,,t) 

MO  =  /  M*(Mh)  \  J1  AM,, 

J  d  Mh 

f  .  An{Mh,t ) 

=  /  e*(M,,)Mh  .J'  AMh 


AM,, 


(24) 


Using  the  power  series  approximation  equation  (19)  together  with 
equations  (22)  and  (23),  the  contribution  to  the  cosmic  SFR  density 
from  haloes  of  mass  M/,  (the  integrand  of  equation  (24))  is  given  by 


dp* 

m, 



1  dMtj 

Mi,  d  t  \ 


Mh 


d n(Mh,  t ) 
dMh 



46230.9p0 

MhStVhT 


x  exp 


232382 
S  p/3,8/3 


(25) 


The  cosmological  constant  term  enters  through  both  the  multiplier 
and  exponential  terms,  with  a  balance  that  depends  on  the  halo  mass 
through  S  (see  equation  17).  While  smaller  haloes  are  more  abundant 
than  large  objects,  a  smaller  fraction  of  the  inflowing  material  is 
converted  into  stars.  As  a  result,  the  SFR  density  is  dominated 
by  the  largest  haloes  in  which  star  formation  is  able  to  proceed 
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Figure  5.  The  predicted  SFR  history  of  the  Universe,  and  the  expected 
influence  of  the  cosmological  constant  using  the  simple  model  developed  in 
Section  2.5.  Coloured  lines  show  the  contributions  from  dark  matter  haloes 
of  different  masses  (per  dex),  using  the  star  formation  efficiency  described  by 
equation  (25).  The  total  SFR  for  the  ACDM  universe  calculated  numerically 
is  shown  in  blue.  An  Einstein-de  Sitter  universe  is  shown  in  orange.  The 
integrated  SFR  calculated  using  the  approximation  of  equations  (24)  and 
(25)  is  shown  with  a  dashed  green  line.  The  bottom  panel  shows  the  ratio 
at  a  given  time.  The  predicted  suppression  of  SFR  due  to  A  at  the  present 
time  is  19  per  cent.  At  t  30  Gys  the  predicted  SFR  density  for  the  EdS 
model  is  double  than  ACDM,  and  ~6  times  higher  at  t  =  50  Gyr.  The 
approximation  of  equations  (24)  and  (25)  ceases  to  work  for  t>  25  Gyr. 

without  generating  efficient  BH  feedback.  The  smaller  haloes  only 
contribute  significantly  at  very  early  times,  when  the  abundance  of 
larger  objects  is  strongly  suppressed  by  the  exponential  term.  We 
see  therefore  that  the  level  of  suppression  expected  for  r«1012  Mg 
haloes  is  representative  of  most  of  the  SFR  in  the  Universe. 

The  predictions  for  the  contributions  of  different  halo  masses 
are  shown  in  Fig.  5,  together  with  the  total  expected  cosmic  SFR 
density,  for  the  two  cosmologies  that  we  consider  in  this  paper.  We 
will  compare  this  approximation  in  Section  5  to  the  results  from  the 
eagle  simulations. 

3  THE  EAGLE  SIMULATIONS 

The  simple  analytic  model  provides  a  basis  for  interpreting  the 
results,  but  it  is  highly  simplified.  We  therefore  compare  the  analytic 
model  to  numerical  hydrodynamic  simulations  based  on  the  eagle 
project.  The  eagle  simulation  suite3  (Crain  et  al.  2015;  Schaye  et  al. 
2015)  consists  of  a  large  number  of  cosmological  hydrodynamical 
simulations  that  include  different  resolutions,  simulated  volumes, 
and  physical  models.  These  simulations  use  advanced  smoothed 

3http://www.eaglesim.org 


particle  hydrodynamics  (SPH)  and  state-of-the-art  subgrid  models 
to  capture  the  unresolved  physics.  The  simulation  suite  was  run 
with  a  modified  version  of  the  gadget-3  SPF1  code  (last  described 
by  Springel  2005)  and  includes  a  full  treatment  of  gravity  and 
hydrodynamics.  The  calibration  strategy  is  described  in  detail  by 
Crain  et  al.  (2015),  who  also  presented  additional  simulations  to 
demonstrate  the  effect  of  parameter  variations. 

The  halo  and  galaxy  catalogues  for  more  than  105  simulated 
galaxies  of  the  main  eagle  simulations  with  integrated  quantities 
describing  the  galaxies,  such  as  stellar  mass,  SFRs.  metallicities,  and 
luminosities,  are  available  in  the  eagle  data  base4  (McAlpine  et  al. 
2016).  A  complete  description  of  the  code  and  physical  parameters 
used  can  be  found  in  Schaye  et  al.  (2015). 

The  eagle  reference  simulations  used  cosmological  parameters 
measured  by  Planck  Collaboration  XVI  (2014).  In  this  paper  we 
introduce  three  main  eagle  simulations  that  use  the  same  calibrated 
sub-grid  parameters  as  the  reference  model,  but  change  the  cosmo¬ 
logical  model  by  setting  the  cosmological  constant  to  zero,  and/or 
removing  feedback  from  BHs.  The  values  of  the  cosmological  pa¬ 
rameters  used  for  the  simulations  are  listed  in  Table  1 .  The  values 
of  other  relevant  parameters  adopted  by  all  simulations  featured  in 
this  study  are  listed  in  Table  2.  Together,  these  parameters  deter¬ 
mine  the  dynamic  range  and  resolution  that  can  be  achieved  by  the 
simulations. 

Fig.  6  shows  the  projected  gas  density  for  the  ACDM  and  EdS 
cosmological  models  both  at  the  present  day  and  into  the  future. 
At  t  =  13.8  Gyr,  the  general  appearance  of  both  models  is  similar, 
but  over  the  next  6.8  Gyr,  the  effect  of  A  becomes  more  significant 
slowing  down  the  growth  of  structure. 

3.1  Subgrid  models 

Processes  that  are  not  resolved  by  the  simulations  are  implemented 
as  subgrid  physical  models;  they  depend  solely  on  local  interstellar 
medium  (ISM)  properties.  A  full  description  of  these  subgrid  models 
can  be  found  in  Schaye  et  al.  (2015).  In  summary: 

(i)  Radiative  cooling  and  photoheating  are  implemented  element 
by  element  as  in  Wiersma,  Schaye  &  Smith  (2009a),  including 
the  11  elements  found  to  be  important,  namely,  H,  He,  C,  N,  O, 
Ne,  Mg,  Si,  S,  Ca,  and  Fe.  Hydrogen  reionization  is  implemented 
by  switching  on  the  full  Haardt  &  Madau  (2001)  background  at 
the  proper  time  corresponding  to  redshift  z  =  11.5  in  our  ACDM 
Universe. 

(ii)  Star  formation  is  implemented  stochastically  following  the 
pressure-dependent  Kennicutt-Schmidt  relation  as  in  Schaye  & 
Dalla  Vecchia  (2008).  Above  a  metallicity-dependent  density 
threshold  nalAJu(Z),  which  is  designed  to  track  the  transition  from 
a  warm  atomic  to  an  unresolved,  cold  molecular  gas  phase  (Schaye 
2004),  gas  particles  have  a  probability  of  forming  stars  determined 
by  their  pressure. 

(iii)  Time-dependent  stellar  mass  loss  due  to  winds  from  mas¬ 
sive  stars  and  AGB  stars,  core  collapse  supernovae  and  type  la 
supernovae,  is  tracked  following  Wiersma  et  al.  (2009b). 

(iv)  Stellar  feedback  is  treated  stochastically,  using  the  thermal 
injection  method  described  in  Dalla  Vecchia  &  Schaye  (2012). 

(v)  Seed  BHs  of  mass  M  =  1.48  x  105  Mg,  are  placed  in  haloes 
with  a  mass  greater  than  1.48  x  1010Mq  and  tracked  following 
the  methodology  of  Springel,  Di  Matteo  &  Hernquist  (2005)  and 


4http://www.eaglesim.org/database.php 
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Table  2.  Box-size,  number  of  particles,  initial  baryonic,  and  dark  matter  particle  mass,  comoving  and  Plummer-equivalent  gravitational  softening,  inclusion 
of  AGN  feedback,  cosmological  model,  and  Hubble  parameter  for  the  eagle  simulations  used  in  this  paper.  Values  in  bold  show  differences  with  respect  to  the 
Ref  simulation.  The  three  bottom  small  box  models  were  used  for  convergence  tests. 


Identifier 

L 

(cMpc) 

N 

mgas 

(M©) 

mDM 

(M0) 

^com 

(ckpc) 

£prop 

(pkpc) 

AGN 

Cosmology 

h 

ACDM  (Ref) 

50 

2  x 

7523 

1.81 

X 

10s 

9.70  x 

106 

2.66 

0.70 

Yes 

Planck  14 

0.6777 

ACDM  (No 
AGN) 

50 

2  x 

7523 

1.81 

X 

106 

9.70  x 

106 

2.66 

0.70 

No 

Planck  14 

0.6777 

EdS 

50 

2  x 

7523 

1.81 

X 

106 

9.70  x 

106 

2.66 

0.70 

Yes 

EdS 

0.3754 

EdS  (No 

AGN) 

50 

2  x 

7523 

1.81 

X 

106 

9.70  x 

106 

2.66 

0.70 

No 

EdS 

0.3754 

A  =  0  L12 

ML3754 

12.50 

2  x 

1883 

1.81 

X 

106 

9.70  x 

106 

2.66 

0.70 

Yes 

EdS 

0.3754 

A  =  0  L12 

h0_6777 

8.43 

2  x 

1883 

1.81 

X 

106 

9.70  x 

106 

1.79 

0.70 

Yes 

EdS 

0.6777 

A  =  0  L12 

h0_4716 

10.73 

2  x 

1883 

1.81 

X 

106 

9.70  x 

106 

2.28 

0.70 

Yes 

EdS 

0.4716 

Booth  &  Schaye  (2009).  Accretion  on  to  BHs  follows  a  modified 
version  of  the  Bondi-Hoyle  accretion  rate  which  takes  into  account 
the  circularization  and  subsequent  viscous  transport  of  infalling 
material,  limited  by  the  Eddington  rate  as  described  by  Rosas- 
Guevara  et  al.  (2015).5  Additionally,  BHs  can  grow  by  merging 
with  other  BHs  as  described  in  Schaye  et  al.  (2015)  and  Salcido 
et  al.  (2016). 

(vi)  Feedback  from  AGN  is  implemented  following  the  stochas¬ 
tic  heating  scheme  described  by  Schaye  et  al.  (2015).  Similar  to  the 
supernova  feedback,  a  fraction  of  the  accreted  gas  on  to  the  BH  is 
released  as  thermal  energy  with  a  fixed  heating  temperature  into  the 
surrounding  gas  following  Booth  &  Schaye  (2009). 

For  the  eagle  simulations,  the  subgrid  parameters  were  calibrated 
to  reproduce  three  properties  of  galaxies  at  redshift  z  =  0:  the  galaxy 
stellar  mass  function  (GSMF),  the  galaxy  size-stellar  mass  relation, 
and  the  BH  mass-stellar  mass  relation.6  The  calibration  strategy  is 
described  in  detail  by  Crain  et  al.  (2015),  who  explore  the  effect  of 
parameter  variations. 

3.2  Halo  and  galaxy  definition 

Haloes  were  identified  running  the  ’Friends-of-Friends’  (FoF)  halo 
finder  on  the  dark  matter  distribution,  with  a  linking  length  equal  to 
0.2  times  the  mean  inter-particle  spacing.  Galaxies  were  identified 
as  self-bound  over-densities  within  the  FoF  group  using  the  subfind 
algorithm  (Springel  et  al.  2001;  Dolag  et  al.  2009).  A  ‘central’ 
galaxy  is  the  substructure  with  the  largest  mass  within  a  halo.  All 
other  substructures  within  a  halo  are  ‘satellite’  galaxies. 

Comparing  haloes  from  simulations  with  different  cosmologies  is 
not  a  well-defined  task,  as  halo  masses  are  usually  defined  in  terms 
of  quantities  that  depend  on  the  specific  cosmological  parameters. 
Typically,  this  is  done  by  growing  a  sphere  outwards  from  the 
potential  minimum  of  the  dominant  dark  matter  sub-halo  out  to  a 
radius  where  the  mean  interior  density  equals  a  fixed  multiple  of 
the  critical  or  mean  density  of  the  Universe,  causing  an  artificial 
‘pseudo-evolution"  of  dark  matter  haloes  by  changing  the  radius  of 
the  halo  (Diemer,  More  &  Kravtsov  2013).  Star  formation,  however, 
is  governed  by  the  amount  of  gas  that  enters  these  haloes  and  reaches 

5The  eagle  simulation  does  not  include  a  boost  factor  the  accretion  rate  of 
BHs  to  account  for  an  unresolved  clumping  factor. 

6BH  feedback  efficiency  left  unchanged  from  Booth  &  Schaye  (2009). 


their  central  regions.  Wetzel  &  Nagai  (2015)  show  that  the  growth 
of  dark  matter  haloes  is  subject  to  this  ‘pseudo-evolution’,  whereas 
the  accretion  of  gas  is  not.  Because  gas  is  able  to  cool  radiatively,  it 
decouples  from  dark  matter,  tracking  the  accretion  rate  near  a  radius 
of  Rioap,  the  radius  within  which  the  mean  density  is  200  times  the 
mean  density  of  the  universe,  p.  As  we  try  to  connect  the  accretion 
of  dark  matter  haloes  to  star  formation,  we  define  halo  masses  as 
the  total  mass  within  Rimp, 

M200?  =  200  ^4^.  (26) 

Additionally,  as  p  =  is  given  in  comoving  coordinates, 

the  mean  density  of  the  universe  remains  constant  in  time  for  each 
cosmological  model. 

Following  Schaye  et  al.  (2015)  and  Furlong  et  al.  (2015),  galaxy 
stellar  masses  are  defined  as  the  stellar  mass  associated  with  the 
subhalo  within  a  3D  30  proper  kilo  parsec  (pkpc)  radius,  centred 
on  the  minimum  of  the  subhalo’s  centre  of  gravitational  potential. 
This  definition  is  equivalent  to  the  total  subhalo  mass  for  low-mass 
objects,  but  excludes  diffuse  mass  around  very  large  subhaloes, 
which  would  contribute  to  the  intracluster  light. 

3.3  Continuing  the  simulations  into  the  future 

As  A  continues  driving  the  accelerated  expansion  of  the  universe, 
the  linear  growth  of  density  perturbations,  £)(/),  is  suppressed  (see 
equation  15).  Further  insight  can  be  obtained  if  we  analyse  the  evo¬ 
lution  of  the  potential  perturbations  given  by  the  perturbed  Poisson 
equation  for  an  expanding  space, 

V2cj>  =  4nGppa2D80,  (27) 

where  the  Laplace  operator  is  with  respect  to  comoving  coordi¬ 
nates,  and  the  mean  density  pp  is  given  in  proper  coordinates.  As 
Pp  evolves  oc  or3,  it  follows  that  V2|t>  oc  D/a.  Using  equations  (10) 
and  (15),  we  can  see  that  for  an  EdS  universe,  both  D  and  a  are 
oc  t213  and  the  potentials  are  expected  to  stop  evolving  (they  are 
frozen  in).  On  the  other  hand,  the  suppression  of  growth  of  density 
perturbations  due  to  a  cosmological  constant  causes  a  decay  in  the 
potentials  as  the  universe  expands.  As  shown  in  Fig.  2,  according  to 
linear  theory,  these  two  scenarios  have  comparable  growth  factors 
at  the  present  time  (^lOpercent  difference;  see  equation  15),  but 
the  difference  becomes  increasingly  important  in  the  future.  Fur¬ 
thermore,  star  formation  is  expected  to  eventually  exhaust  the  finite 
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Figure  6.  The  evolution  of  the  projected  gas  density  for  each  eagle  model  centred  on  the  most  massive  halo  at  the  present  time  ( t  =  13.8  Gyr).  The  length 
of  each  image  is  43  (proper)  Mpc  on  a  side,  to  highlight  the  difference  on  cosmic  expansion.  Left:  ACDM  universe.  Right:  EdS  universe.  Top:  Cosmic  time 
t  =  13.8  Gyr.  Bottom:  Cosmic  time  t  =  20.7  Gyr.  The  colour  coding  represents  the  (proper)  surface  gas  density  projected  along  the  line  of  sight.  At  t  =  13.8  Gyr, 
the  general  appearance  of  both  models  is  similar,  as  the  phases  of  the  initial  fluctuations  are  the  same.  Over  the  next  6.8  Gyr,  the  effect  of  A  becomes  more 
significant,  slowing  down  the  growth  of  structure  compared  to  the  EdS  model. 


reservoir  of  cold  gas  in  galaxies,  shutting  off  the  production  of  stars 
in  the  universe  forever  (e.g.  Fukugita,  Hogan  &  Peebles  1998;  Loeb, 
Batista  &  Sloan  2016). 

In  order  to  study  the  impact  of  A  in  galaxy  formation  beyond  the 
present  day,  and  hence  explore  the  uniqueness  of  the  present  epoch, 
and  in  order  to  determine  the  total  mass  of  stars  ever  produced  by  the 
universe,  we  allow  the  simulations  to  run  into  the  future,  i.e.  t  >  t0 
(e.g.  Barnes  et  al.  2005;  Loeb  et  al.  2016).  The  subgrid  models  for 
star  formation,  stellar  mass  loss,  stellar  feedback,  BH  seeding,  and 
feedback  from  AGN  were  kept  as  described  in  Section  3.1  as  the 
simulations  ran  into  the  future.  On  the  other  hand,  as  there  is  no 
information  about  the  UV  and  X-ray  background  radiation  from 
quasars  and  galaxies  into  the  future,  for  simplicity,  we  assumed  that 
the  background  radiation  freezes  out,  i.e.  we  kept  its  value  at  t  —  to 
constant  into  the  future.  We  consider  this  to  be  a  good  simplification 


as  the  UV  background  only  affects  star  formation  in  very  low  mass 
haloes,  and  hence  does  not  affect  the  cosmic  SFR  at  late  times  (e.g. 
Schaye  et  al.  2010). 


4  SIMULATIONS  RE-SCALING 

In  this  section  we  describe  our  simulation  re-scaling  strategy.  At 
early  epochs,  the  universe  was  matter  dominated,  and  so  we  can 
neglect  the  contribution  of  A.  Hence,  any  universe  with  non-zero 
matter  density,  i.e.  pm>  o  ^  0,  will  be  close  to  an  EdS  universe  at 
early  epochs.  Therefore,  we  can  assume  identical  initial  conditions 
for  all  cosmological  models  of  interest  here. 

The  initial  conditions  for  the  reference  ACDM  model  were  cre¬ 
ated  in  three  steps.  First,  a  particle  load,  representing  an  unperturbed 
homogeneous  periodic  universe,  was  produced.  Secondly,  a  realiza- 
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Table  3.  Parameters  re-scaled  in  the  initial  conditions.  Hat  notation  indi¬ 
cates  parameters  for  our  Universe. 


Parameter 

Units 

Re-scaling  factor 

Box  size 

cMpc/U1 

(h~lh)  x  (aiflf1) 

Particle  masses 

M0r‘ 

(h-lh) 

Particle  coordinates 

cMpc/U1 

(h~lh)  x  (aiflf1) 

Particle  velocities 

cMpcs^1 

(£na-W!2 

Table  4.  Additional  parameters  re-scaled  in 
indicates  parameters  for  our  Universe. 

the  simulations.  Hat  notation 

Parameter 

Units 

Re- scaling  factor 

Comoving  softening 

ckpc  h~ 1 

(h~xh)  x  (uia]"1) 

Max  softening 

pkpc/i-1 

(/A1/;) 

Seed  BH  mass 

Mq/U1 

(/A1/!) 

Min  Mfof  for  new  BH 

Mq/U1 

(/A'/i) 

ACDM  Redshift  z 

109  876543  2  1  0 


CO 

I 


<C3 


* 

•a 


CO 


o 


o 


Age  of  the  universe  [Gyr] 


tion  of  a  Gaussian  random  density  field  with  the  appropriate  linear 
power  spectrum  was  created  over  the  periodic  volume.  Thirdly, 
the  displacements  and  velocities,  consistent  with  the  pure  growing 
mode  of  gravitational  instability,  were  calculated  from  the  Gaus¬ 
sian  realization  and  applied  to  the  particle  load  producing  the  initial 
conditions.  The  initial  density  perturbation  power  spectrum  is  com¬ 
monly  assumed  to  be  a  power  law,  i.e.  Pt  ( k )  oc  kn‘ .  From  the  Planck 
results  (Planck  Collaboration  XVI  2014),  the  spectral  index  ns,  has 
a  value  of  ns  =  0.9611.  A  transfer  function  with  the  cosmological 
parameters  shown  in  Table  1  was  generated  using  CAMB  (version 
Jan_12;  Lewis,  Challinor  &  Lasenby  2000).  The  linear  matter  power 
spectrum  was  generated  by  multiplying  the  initial  power  spectrum 
by  the  square  of  the  dark  matter  transfer  function  evaluated  at  the 
present  day  t  =  t0,  i.e.  P(k,  t)  =  Pi(k)T2(k)D2(t)3 

The  eagle  version  of  gadget  uses  an  internal  system  of  units  that 
includes  both  comoving  coordinates  and  the  dimensionless  Hubble 
parameter,  h.  For  the  alternative  cosmological  models,  we  have  the 
freedom  to  choose  the  present  time  t0  for  each  simulation,  and  we  re¬ 
scale  all  the  initial  condition  such  that  they  are  identical  in  physical 
’h-free'  units  at  an  early  time  t\  =  11.98Myr.  Table  3  shows  the 
parameters  that  have  been  re-scaled  in  the  initial  conditions. 

The  same  tables  of  radiative  cooling  and  photoheating  rates  as  a 
function  of  density  and  temperature  were  used  for  all  cosmological 
models.  The  corresponding  redshifts  for  the  cooling  tables  were  re¬ 
scaled  such  that  they  correspond  to  the  same  cosmic  time  for  each 
cosmology.  That  is,  using  equation  (8),  we  find  the  scale  factor,  a, 
for  which  the  alternative  cosmology  satisfies 

t{a)  =  t(a).  (28) 


The  average  baryonic  density  £2b  has  been  re-scaled  in  such  a  way 
that  the  baryon  fraction  (fb  —  I2b/I2m)  is  equal  in  both  cosmologies, 
i.e. 


G  |, 


(29) 


Table  4  shows  additional  parameters  that  have  been  re-scaled  to  be 
equivalent  in  /?-free  physical  units.  Finally,  hydrogen  and  He  ii  reion¬ 
ization  were  also  re-scaled  in  such  a  way  that  redshifts  correspond 
to  the  same  cosmic  time. 


7The  CAMB  input  parameter  file  and  the  linear  power  spectrum  are  available 
at  http://eagle.strw.leidenuniv.nl/ 


Figure  7.  Global  SFR  density  for  three  EdS  models  scaled  by  the  ratio  of 
the  initial  scale  factors  for  each  model.  The  initial  conditions  for  each  model 
have  been  re-scaled  such  that  the  time  at  which  we  start  the  simulations 
remains  unchanged,  i.e.  t(a\)  =  11.98  Myr. 

In  order  to  demonstrate  that  this  re-scaling  strategy  works  cor¬ 
rectly,  Fig.  7  shows  the  global  SFR  density  for  the  three  small  box 
EdS  simulations  used  for  convergence  (see  Table  2).  They  each 
represent  the  same  physical  scenario,  but  choose  a  different  proper 
time  to  be  ’today’,  t0.  This  has  the  effect  of  altering  the  values  of 
the  Hubble  parameter  h  and  the  redshift  of  the  initial  conditions, 
so  that  the  simulations  begin  at  proper  time  t(ai)  =  11.98 Myr  in 
all  models.  Despite  the  small  size  of  the  simulation  boxes  (hence 
the  noisy  curves),  the  figure  shows  consistent  SFRs  as  a  function  of 
cosmic  time  for  the  three  models.  Therefore,  our  re-scaling  strategy 
allows  us  to  simulate  any  cosmological  model,  regardless  of  the 
value  of  h. 


5  RESULTS:  THE  EVOLUTION  OF  STAR 
FORMATION 

5.1  The  past  history  of  the  cosmic  SFR 

Fig.  8  shows  the  global  SFR  density  as  a  function  of  cosmic  time  for 
our  simulation  models.  For  comparison,  observations  from  Cucciati 
et  al.  (2012)  [FUV],  Bouwens  et  al.  (2012)  [UV],  Robertson  et  al. 
(2013)  [UV],  and  Burgarella  et  al.  (2013)  [FUV  +  FIR]  are  shown  as 
well.  Solid  lines  in  the  figure  show  the  evolution  of  the  (comoving) 
cosmic  SFR  density  for  the  reference  ACDM  eagle  run  (blue)  and 
for  an  EdS  universe  (orange).  Dashed  lines  show  simulation  models 
without  feedback  from  AGN.  Dotted  lines  show  the  prediction  for 
the  cosmic  SFR  density  using  equation  (25).  We  focus  first  on 
the  evolution  of  the  models  up  to  the  present  age  of  the  universe, 
t  =  t0  =  13.8  Gyr. 

In  linear  time,  the  SFR  rises  very  rapidly  and  most  of  the  plot  is 
dominated  by  the  slow  decline  (for  an  example,  see  Furlong  et  al. 
2015).  Hence,  in  order  to  emphasize  the  growth  and  decline  of 
the  SFR,  and  to  reproduce  the  familiar  shape  of  the  star  formation 
history  (Madau  &  Dickinson  2014),  the  horizontal  time  axis  has 
been  plotted  in  a  logarithmic  scale  for  t  <  8  Gyr.  In  order  to  explore 
the  SFR  in  detail  at  the  present  epoch  and  into  the  future,  the 
horizontal  time  axis  changes  to  a  linear  scale  for  t  >  8  Gyr.  The 
black  vertical  dotted  line  shows  the  transition  from  logarithmic  to 
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Figure  8.  Global  SFR  densities.  The  eagle  reference  and  the  EdS  models  are  shown  in  solid  blue  and  orange  lines,  respectively.  Observational  data  from 
Cucciati  et  al.  (2012)  [FUV],  Bouwens  et  al.  (2012)  [UV],  Robertson  et  al.  (2013)  [UV],  and  Burgarella  et  al.  (2013)  [FUV  +  FIR]  are  shown  as  symbols. 
The  eagle  reference  simulation  (solid  blue)  reproduces  the  shape  of  the  observed  SFR  density  remarkably  well,  with  a  small  offset  of  0.2  dex  at  1  >  2  Gyr 
(for  clarity,  all  eagle  models  have  been  shifted  by  0.2  dex).  The  analytical  model  of  equations  (24)  and  (25)  is  shown  with  dotted  lines  for  both  cosmologies. 
Power-law  fittings  for  the  SFR  density  for  t  >  8  Gyr,  as  per  equation  (30),  for  the  ACDM  and  EdS  models  are  shown  in  pink  and  red  lines,  respectively.  The 
horizontal  time  axis  has  been  plotted  in  a  logarithmic  scale  for  t  <  8  Gyr  changing  to  a  linear  scale  for  l  >  8  Gyr.  The  black  vertical  dotted  line  shows  the 
transition  from  logarithmic  to  linear  scale. 


linear  scale.  For  reference,  the  redshifts,  z,  for  an  observer  at  tg  in  the 
ACDM  universe,  are  given  along  the  top  axis.  As  discussed  in  detail 
in  Furlong  et  al.  (2015),  the  reference  simulation  (solid  blue  line) 
reproduces  the  shape  of  the  observed  SFR  density  remarkably  well, 
with  a  small  offset  of  0.2  dex  at  t  >  2  Gyr.  While  the  simulations 
agree  reasonably  well  with  the  observational  data  at  redshifts  above 
3,  we  caution  that  these  measurements  are  more  uncertain. 

Remarkably,  the  shape  of  the  cosmic  SFR  history  is  very  sim¬ 
ilar  for  both  the  ACDM  and  EdS  models:  the  SFR  density  peaks 
^3.5  Gyr  after  the  Big  Bang  and  declines  slowly  thereafter.  The 
similarity  of  the  universes  prior  to  the  peak  is  expected,  since  the 
A  term  in  the  Friedman  equation  is  sub-dominant  in  both  cases.  At 
later  times,  however,  we  might  naively  have  expected  the  decline  to 
be  more  pronounced  in  the  ACDM  cosmology,  since  the  growing 
importance  of  the  A  term  slows  the  growth  of  density  perturbations. 

From  Fig.  2,  the  linear  growth  factors  of  the  two  cosmological 
models  differ  by  sa  10  per  cent  at  the  present  time,  and  so  we  might 
have  expected  a  similar  difference  in  the  (comoving)  cosmic  SFR 
density  («15  percent,  read  from  Fig.  8).  This  naive  expectation  is 
not  borne  out  because  of  the  complexity  of  the  baryonic  physics. 
Because  of  stellar  and  AGN  feedback,  haloes  have  an  ample  reser¬ 
voir  of  cooling  gas  that  is  able  to  power  further  star  formation, 
regardless  of  the  change  in  the  cosmic  halo  growth  rate. 

Our  simulation  demonstrates  that  the  existence  of  A  does  play 
a  small  role  in  determining  the  (comoving)  cosmic  SFR  density. 
Flowever,  these  differences  are  minor.  In  order  to  put  the  differences 


into  context,  we  compare  with  a  pair  of  simulations  in  which  the  BH 
feedback  is  absent.  These  runs  are  shown  as  dashed  lines  in  the  plot. 
We  focus  here  on  the  behaviour  before  t  =  13.8  yr.  As  can  be  seen, 
the  absence  of  AGN  feedback  has  a  dramatic  effect  on  the  shape 
of  the  cosmic  star  formation  density  (Schaye  et  al.  2010;  van  de 
Voort  et  al.  2011).  Interestingly,  however,  while  the  normalization 
of  the  SFR  density  is  considerably  higher,  the  time  of  the  peak  is 
similar.  BH  feedback  is  not  solely  responsible  for  the  decline  in  star 
formation  after  t  3.5  Gyr.  This  hints  that  the  existence  of  the  peak 
results  from  the  interaction  of  the  slowing  growth  rates  of  haloes 
(after  the  peak)  and  the  star  formation  time-scale  (set  by  the  ISM 
physics)  which  limits  the  rate  at  which  the  galaxy  can  respond  to 
convert  in-falling  material  into  stars  (before  the  peak). 

The  SFR  history  predicted  by  the  simple  model  developed  in 
Section  2.5  (dotted  curves)  is  in  remarkable  agreement  with  the 
observational  data  and  predicts  the  relative  difference  of  the  cos¬ 
mological  simulations,  both  at  the  present  time  and  into  the  future. 
We  want  to  emphasize  that  the  model  is  not  a  parametric  fit  to  the 
data,  but  rather  an  analytical  model  derived  from  a  simple  relation 
of  star  formation  to  halo  mass  accretion. 


5.2  The  future  of  the  cosmic  SFR  ( t  >  13.8  Gyr) 

In  order  to  explore  whether  the  relative  SFR  densities  will  diverge  as 
the  impact  of  A  becomes  more  pronounced,  we  ran  the  simulations 
for  both  cosmological  models  into  the  future  (i.e.  beyond  a  cosmic 
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Table  5.  Power-law  parameter  fitting  for  the  median  SFR  shown  in  Fig.  8. 


Model 

a 

c 

k 

(Mg  yr_1  cMpc-3) 

(Gyr-1) 

ACDM 

5.42 

0.81 

2.77 

EdS 

3.99 

0.94 

2.41 

time  of  t  =  13.8  Gyr).  As  the  simulations  run  into  the  future,  the 
small  differences  seen  at  t  =  13.8  Gyr  become  larger,  reading  an 
«40  per  cent  difference  at  t  =  1.5  x  to  =  20.7  Gyr,  which  is  in 
agreement  with  the  predictions  shown  in  Fig.  5. 

The  decline  in  the  SFR  density  can  be  approximated  by  a  power 
law  for  both  the  reference  ACDM  and  EdS  models  (red  and  pink 
dashed  lines), 

A  =  a{ctyk,  (30) 

with  the  parameters  a,  c,  and  k  given  in  Table  5.  We  used  the  reduced 
chi-square  statistic  for  goodness-of-fit  testing.  We  use  this  fitting 
function  to  extrapolate  the  results  from  the  simulations  further  into 
the  future. 

We  note  here  a  striking  feature  of  the  universes  with  no  BF1 
feedback:  the  SFR  increases  again  in  the  future,  for  both  the  ACDM 
and  EdS  cosmologies.  In  Section  5.4  we  discuss  how  this  effect 
originates  from  massive  galaxies  (M*  >  1011  Mg)  that  rejuvenate 
in  the  future.  As  it  is  not  clear  from  the  simulation  whether  the  SFR 
will  continue  to  rise,  or  when  it  would  start  declining,  we  have  not 
fitted  any  functional  form  to  the  "No  AGN'  models. 

5.3  The  stellar  mass  density 

To  study  the  build-up  of  stellar  mass,  we  present  the  growth  in 
stellar  mass  density,  pt,  across  cosmic  time  in  Fig.  9.  The  colour 
coding  is  the  same  as  in  Fig.  8.  The  lower  panel  shows  the  ratio 
of  the  stellar  density  compared  to  the  reference  ACDM  model. 
Furlong  et  al.  (2015)  show  that  the  reference  eagle  simulation  is 
in  good  agreement  with  the  observed  growth  of  stellar  mass  across 
cosmic  time.  In  contrast  to  Furlong  et  al.  (2015),  where  the  stellar 
mass  density  was  obtained  from  aperture  measurements  to  facilitate 
comparison  with  observations,  we  calculate  pt  by  integrating  the 
SFR  density  from  Fig.  8, 

P*  =  [  P,dt',  (31) 

Jo 

in  order  to  provide  an  estimate  of  the  total  mass  of  stars  produced  by 
the  universe.8  For  the  models  with  AGN  feedback,  we  have  extrap¬ 
olated  the  power-law  fit  described  in  Section  5.1  far  into  the  future, 
up  to  10  trillion  years,  and  considered  this  as  the  ’Total  stellar  mass 
density’  of  the  universe.  As  suggested  by  the  analytic  model  in 
equation  (25),  in  universes  with  feedback  from  star  formation  and 
AGN,  the  cosmic  SFR  density  is  expected  to  continue  decreasing 
into  the  future.  At  late  times,  the  SFR  becomes  orders  of  magnitude 
lower  than  that  of  the  peak  at  3.5  Gyr.  Fig.  9  shows  that  the  total 
stellar  mass  density  is  dominated  by  the  contribution  from  the  peak 
in  star  formation  and  reaches  a  plateau  at  t  m  to-  Hence,  the  formal 
uncertainties  in  the  extrapolation  into  the  future  are  unimportant 
for  the  predicted  total  stellar  mass  density.  The  right-hand  axis  of 
Fig.  9  represents  the  percentage  of  the  total  stellar  density  com¬ 
pared  to  ACDM.  For  the  reference  ACDM  model,  the  universe  has 

8  As  we  are  interested  in  total  mass  of  stars  produced  by  the  universe,  equation 
(31)  ignores  stellar  mass  loss. 


ACDM  Redshift  z 

10  8765432  1  0 


Figure  9.  The  stellar  mass  density  as  a  function  of  time  in  the  eagle  simula¬ 
tion  models.  The  colour  coding  is  the  same  as  in  Fig.  8.  The  right-hand  axis 
represents  the  percentage  of  the  total  stellar  density  compared  to  the  ACDM 
model.  Both  the  ACDM  and  EdS  models  have  already  produced  most  of  the 
stars  in  the  universe  by  the  present  day,  building  up  very  little  stellar  mass 
into  the  future.  The  models  without  AGN  feedback  quickly  deviate  from  the 
reference  model,  producing  almost  twice  the  mass  in  stars  by  the  end  of  the 
simulation  (~20.7  Gyr).  The  figure  shows  that  the  effect  of  dark  energy  on 
the  overall  star  formation  is  negligible. 

already  produced  most  (^88 percent)  of  its  eventual  stellar  mass 
by  the  present  day,  adding  up  very  little  stellar  mass  into  the  future. 
Although  there  is  no  A  to  slow  down  the  formation  of  cosmic  struc¬ 
ture,  the  EdS  model  closely  resembles  a  ACDM  cosmology  with 
only  ^15  percent  more  stellar  mass  produced. 

As  discussed  in  Section  3.3,  a  universe  with  A  has  a  very  differ¬ 
ent  expansion  history  compared  to  one  without  dark  energy.  This 
produces  a  different  growth  of  density  perturbations,  in  particular 
into  the  future  (see  Figs  1  and  2).  Nevertheless,  as  seen  in  Fig.  9, 
since  both  cosmologies  have  already  produced  most  of  the  stars 
in  the  universe  by  the  present  day,  when  the  contribution  of  A  is 
becoming  increasingly  important,  the  effect  of  dark  energy  on  the 
overall  star  formation  is  negligible. 

In  contrast,  the  models  without  feedback  from  BHs  (dashed 
curves)  quickly  deviate  from  the  reference  model,  starting  from 
t  ~  1  Gyr,  producing  almost  twice  the  mass  in  stars  by  the  end  of 
the  simulations,  20.7  Gyr  after  the  Big  Bang. 

5.4  Other  galaxy  population  properties 

In  this  section  we  will  compare  the  galaxy  population  properties  of 
the  two  simulation  models  at  the  present  time  and  into  the  future. 
In  particular,  we  compare  the  GSMF  and  the  specific  star  formation 
rate  (SSFR)  of  galaxies.  To  compare  with  observational  data,  in  each 
of  the  figures  discussed  below,  the  left-hand  panel  shows  properties 
at  /  =  12.5  Gyr,  equivalent  to  redshift  2  =  0.1  for  an  observer  at 
the  present  time  in  a  ACDM  universe.  The  right-hand  panel  shows 
the  same  property  but  at  t  =  1 .5  x  to  =  20.7  Gyr.  To  guide  the  eye, 
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Figure  10.  The  GSMF  at  /  =  12.5  Gyr,  equivalent  to  redshift  z.  =  0.1  for  an  observer  at  the  present  time  in  a  ACDM  universe  (left),  and  t  =  1.5  x  to  (right) 
for  the  eagle  simulation  models.  The  colour  coding  is  the  same  as  in  Fig.  8.  Observational  data  from  Baldry  et  al.  (2012)  and  Li  &  White  (2009)  are  shown 
as  symbols.  The  reference  and  ’No  AGN'  ACDM  models  at  t  =  12.5  Gyr  are  plotted  in  the  right-hand  panel  for  reference  (solid  and  dashed  black  lines, 
respectively).  The  effect  of  dark  energy  on  the  GSMF  is  negligible,  with  very  little  evolution  into  the  future.  As  expected,  the  models  without  AGN  feedback 
predict  a  higher  number  density  of  massive  galaxies  (M*  >  1 01 1  Mq).  This  effect  becomes  more  significant  into  the  future,  going  from  ~0.7  to  swl.l  dex. 


each  property  for  the  reference  ACDM  model  at  t  =  12.5  Gyr  is 
plotted  with  a  black  line  in  each  plot  at  t  =  1.5  x  to  =  20.7  Gyr. 
Finally,  we  have  also  included  the  ratios  of  each  quantity  to  the 
reference  ACDM  model  at  the  corresponding  time  at  the  bottom  of 
each  panel. 


5.4. 1  The  galaxy  stellar  mass  function 

The  effect  of  A  on  the  GSMF  can  be  seen  in  Fig.  10.  The  colour 
coding  is  the  same  as  in  Fig.  8.  For  comparison,  observational  data 
from  Baldry  et  al.  (2012)  and  Li  &  White  (2009)  are  shown  as  well. 
As  discussed  in  Schaye  et  al.  (2015),  the  observed  GSMF  at  redshift 
0.1  was  used  to  infer  the  free  parameters  of  the  subgrid  physics 
used  in  the  simulation.  The  reference  eagle  model  reproduces  the 
shape  of  the  observed  GSMF  reasonably  well,  with  a  slight  under¬ 
abundance  of  galaxies  at  its  knee. 

Fig.  10  shows  that  the  effect  of  A  on  the  GSMF  is  negligi¬ 
ble,  with  very  little  evolution  into  the  future.  The  models  without 
AGN  feedback  predict  a  higher  number  density  of  massive  galaxies 
(M*  >  1011  Mq)  compared  to  the  models  with  AGN  feedback.  This 
effect  becomes  more  significant  into  the  future,  going  from  0.7  to 
1 .5  dex.  The  origin  of  this  difference  is  explored  in  the  next  section. 

5.4.2  Specific  SFRs 

Galaxies  can  be  broadly  classified  into  largely  distinct  star-forming 
and  passive  populations  according  to  their  SSFR, 

M* 

SSFR  =  — (32) 


For  star-forming  galaxies,  there  is  a  well-defined  star-forming  se¬ 
quence,  with  SSFR  observed  to  be  approximately  constant  as  a 
function  of  stellar  mass  (e.g.  Noeske  et  al.  2007;  Karim  et  al.  2011). 
Fig.  11  shows  the  SSFR  for  star-forming  galaxies  in  the  simula¬ 
tions  as  a  function  of  galaxy  stellar  mass  at  the  present  day  and 
into  the  future.  The  colour  coding  is  the  same  as  in  Fig.  8.  The 
horizontal  dotted  lines  correspond  to  the  SSFR  cut  ( 10  2  Gyr-1) 
used  to  separate  star  forming  from  passive  galaxies  in  our  Universe. 
For  comparison,  observational  data  from  Gilbank  et  al.  (2010)  are 
shown  as  well.  Furlong  et  al.  (2015)  show  that  the  SFFR  in  the 
reference  simulations  at  the  present  day  is  similar  to  observations 
in  the  local  universe,  with  an  offset  of  0.3  dex.  This  is  possibly 
consistent  with  the  systematic  uncertainties  in  the  calibration  of  the 
observation  diagnostics.  At  low  masses  there  is  an  increase  in  SSFR 
with  stellar  mass;  however,  this  has  been  found  to  be  a  resolution- 
dependent  effect.  Flence,  we  have  plotted  the  results  with  lighter 
coloured  lines  (similar  to  Schaye  et  al.  2015).  The  models  without 
feedback  from  AGN  have  higher  SSFR  for  M„  >  1010  Mq,  whereas 
the  effect  of  A  on  the  SSFR  of  galaxies  is  negligible. 

Fig.  11  shows  the  galaxy  population  property  that  has  the 
strongest  evolution  into  the  future.  We  find  that  over  the  next  6.8  Gyr 
the  SSFR  will  drop  by  »*0.4  dex. 

Interestingly,  the  models  without  AGN  feedback  predict  an  in¬ 
crease  of  the  SSFR  of  galaxies  with  M*  >  1011  Mq  in  the  future. 
The  figure  shows  that  the  increase  in  SFR  shown  in  Fig.  8,  and 
the  higher  number  density  of  massive  galaxies  in  Fig.  10,  originate 
from  massive  galaxies  that  rejuvenate  in  the  future.  A  plausible 
explanation  for  this  phenomenon  is  that,  according  to  simple  radia¬ 
tive  cooling  models,  in  massive  galaxies  and  clusters,  a  hot  gaseous 
atmosphere  should  lose  energy  by  the  emission  of  radiation,  and 
if  there  is  no  heating  mechanism  to  compensate  the  cooling  (e.g. 
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Figure  11.  The  SSFR  of  star-forming  galaxies  at  t  =  12.5  Gyr,  equivalent  to  redshift  2  =  0.1  for  an  observer  at  the  present  time  in  a  ACDM  universe  (left-hand 
side)  and  (=  1.5  x  to  (right-hand  side)  for  the  eagle  simulation  models.  The  colour  coding  is  the  same  as  in  Fig.  8.  Observational  data  from  Gilbank  et  al. 
(2010)  are  shown  as  symbols.  The  reference  and  'No  AGN’  ACDM  models  at  t  =  12.5  Gyr  are  plotted  in  the  right-hand  panel  for  reference  (solid  and  dashed 
black  lines,  respectively).  The  solid  curves  show  the  median  relation  for  star-forming  galaxies,  defined  as  those  with  an  SSFR  above  the  limit  specified  by  the 
horizontal  dashed  line  (10  2  Gyr  - 1 ).  The  faint  shaded  regions  enclose  the  10th-90th  percentiles  for  the  ACDM  and  EdS  Models.  Lines  are  light  coloured  when 
the  stellar  mass  falls  below  that  corresponding  to  100  baryonic  particles,  to  indicate  that  resolution  effects  will  be  important.  The  figure  shows  that  the  effect 
of  dark  energy  on  the  GSMF  is  negligible.  The  models  without  AGN  feedback  predict  a  higher  SSFR  for  massive  galaxies  (M*  >  1010  Mq).  The  right-hand 
panel  shows  that  the  overall  SSFR  drops  from  t  =  12.5  Gyr  to  t  =  1.5  x  to-  For  the  'No  AGN’  models,  however,  the  SSFR  increases  for  massive  galaxies 
(M*  >  1011  Mq). 


AGN  feedback),  cooling  flows  should  form  (Fabian  1994;  Peter¬ 
son  &  Fabian  2006),  triggering  star  formation.  This  result  will  be 
explored  in  more  detail  in  a  follow-up  paper. 


6  DISCUSSION  AND  CONCLUSIONS 

In  this  paper,  we  explored  the  dependence  of  the  star  formation 
history  of  the  universe  on  the  existence  of  a  cosmological  constant 
and  feedback  from  accreting  BHs.  We  base  our  results  on  the  eagle 
simulation  code,  which  has  been  shown  to  compare  favourably  to 
observational  data,  and  thus,  to  provide  a  good  description  of  the 
formation  of  galaxies  in  our  Universe.  Feedback  from  supermassive 
BHs  has  been  shown  to  be  a  key  ingredient  in  achieving  this  match 
by  suppressing  star  formation  in  massive  haloes  (e.g.  Bower  et  al. 
2006;  Harrison  2017),  while  the  accelerating  expansion  rate  of  the 
Universe  suppresses  the  accretion  rates  of  haloes  at  late  times  (e.g. 
Jenkins  et  al.  1998;  Huterer  et  al.  2015).  Our  study  allows  us  to 
assess  the  relative  importance  of  these  ingredients. 

The  universes  that  we  consider  are  indistinguishable  at  early 
times.  They  share  a  common  epoch  of  equality  and  recombination, 
and  have  equal  amplitudes  and  spectrum  of  density  fluctuations  at 
early  times.  We  take  care  to  compare  the  evolution  of  models  with 
equivalent  starting  points,  and  to  demonstrate  that  the  simulation 
code  correctly  scales  the  different  values  of  the  present-day  expan¬ 
sion  rate  (Hubble  parameter).  When  comparing  the  universes,  it  is 


important  that  we  compare  properties  at  a  fixed  cosmic  time.  Since 
the  processes  of  stellar  (and  biological)  evolution  provide  a  com¬ 
mon  clock,  independent  of  the  large-scale  cosmological  expansion, 
these  provide  an  astrophysically  relevant  comparison. 

We  have  also  developed  an  analytic  model  derived  from  a  simple 
relation  of  star  formation  to  halo  mass  accretion  rate.  Despite  its 
simplicity,  the  model  reproduces  the  overall  shape  of  evolution  of 
the  cosmic  SFR  density.  The  model  and  the  simulations  allow  us 
to  explore  the  effect  of  a  cosmological  constant  term  on  the  cosmic 
SFR  density. 

Our  main  conclusions  are  as  follows: 

(i)  We  find  that  the  existence  of  the  cosmological  constant  has 
little  impact  on  the  star  formation  history  of  the  Universe.  The  SFR 
is  suppressed  by  r«15  percent  at  the  present  time,  and  we  find  that 
the  properties  of  galaxies  are  almost  indistinguishable  in  the  two 
universes. 

(ii)  To  explore  whether  this  is  due  to  the  relatively  recent  dom¬ 
inance  of  the  dark  energy  density  in  our  Universe,  we  continued 
the  simulations  6.8  Gyr  into  the  future.  Even  after  this  time,  the 
comoving  SFR  densities  differ  only  by  r«40  per  cent.  Clearly,  the 
cosmological  constant  has  only  a  marginal  effect  on  the  stellar  con¬ 
tent  of  the  Universe. 

(iii)  Using  the  analytic  model,  we  can  recognize  that  the  existence 
of  the  peak  in  the  SFR  density  results  from  the  interaction  of  the  star 
formation  efficiency  (set  by  the  ISM  physics)  which  limits  the  rate 
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at  which  the  galaxy  can  respond  to  convert  in-falling  material  into 
stars,  the  relative  abundance  of  efficiently  star-forming  haloes  (i.e. 
of  masses  ^1012  Mg),  and  only  at  late  times,  the  slowing  growth 
rates  of  haloes  due  to  the  cosmological  constant. 

(iv)  By  extrapolating  fits  to  the  evolution  of  the  comoving  SFR 
density  into  the  future,  we  show  that,  in  our  Universe,  more  than 
^88 percent  of  the  stars  that  will  ever  be  produced,  have  already 
been  formed  by  the  present  cosmic  time.  In  the  absence  of  dark  en¬ 
ergy,  only  s=s  15  percent  more  stellar  mass  would  have  been  formed 
in  the  same  time.  The  difference  is  small,  bringing  into  question 
whether  the  ‘coincidence  problem"  (the  comparable  energy  densi¬ 
ties  of  matter  and  dark  energy)  can  be  explained  by  an  anthropic 
argument:  the  existence  of  dark  energy  (at  the  observed  value)  has 
negligible  impact  on  the  existence  of  observers  or  the  ability  of  hu¬ 
manity  to  observe  the  cosmos.  In  Barnes  et  al.  (2018),  we  explore 
this  argument  in  more  detail  by  considering  a  wider  range  of  A 
values,  and  determining  the  likelihood  distribution  of  possible  A 
values  conditioning  the  existence  of  observers. 

(v)  In  comparison,  the  existence  of  BHs  has  a  major  impact  on 
the  Universe.  In  the  absence  of  AGN  feedback,  the  comoving  SFR 
density  is  enhanced  by  a  factor  of  2.5  at  the  present  day. 

(vi)  Even  in  a  universe  without  BHs  or  dark  energy,  we  find  that 
the  comoving  SFR  density  peaks  at  3.5  Gyr  (z  &  2  according  to  a 
present-day  observer  in  our  Universe).  The  decline  in  star  formation 
is  however  slower  at  more  recent  times. 

(vii)  For  hypothetical  universes  without  feedback  from  accreting 
BHs,  there  is  a  comeback  of  SFR,  which  increases  again  in  the  fu¬ 
ture.  This  effect  originates  from  massive  galaxies  (M,  >  1011  Mq) 
that  rejuvenate  as  there  is  no  heating  mechanism  to  compensate  the 
cooling,  in  turn,  triggering  star  formation. 
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